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We observe oscillatory decay in the two-point, non-equal 
time, velocity correlation function of homogeneous, isotropic 
turbulence. We found this through a direct numerical simu- 
lation (DNS) of the three dimensional Navier-Stokes (3 — -D 
NS) equation. We give an approximate analytic theory which 
explains this oscillatory behaviour. The wave- number and fre- 
quency dependent effective viscosity turns out to be complex; 
the imaginary part gives rise to the temporal oscillation. We 
find that, at least for the decay at short times, data collapse 
occur among the inertial range velocity wave-vector modes 
with the long time dynamic exponent z — 2/3, but the time 
period of the temporal oscillation is not universal. 

PACS : 47.27.Gs, 47.27.Eq, 05.45.+b, 05.70.Jk 



In homogeneous, isotropic turbulence the main interest 
is to understand the long-range spatio-temporal correla- 
tions exhibited by the velocity field. Towards this end 
one studies the scaling behaviour of the velocity structure 
functions Sp{l, t) = {[vi{0, 0) — Vi{l, t)]P) with respect to I 
and t. Here Vi is the velocity field, I and t are the spatial 
and temporal separations. Kolmogorov had predicted 
through his dimensional analysis argument that the 
equal time structure function Sp{l,0) = [eiy^^. He had 
assumed that in the inertial range i.e., rjd ^ I ^ L [rjd 
and L are, respectively, dissipation and forcing scales) 
Sp{l,Q) is a function of e (the mean energy dissipation 
rate) and I only. But after anomalous scaling properties 
of Sp{l) were discovered through experiments and simu- 
lations 1^, the importance of L (which is also called the 
integral scale) has been recognised. Now it is known that 
Sp{l) ^ I'^p , when is a monotonically increasing, con- 
vex, nonlinear function of p. The negative correction to 
the exponent d(^p = (^p — p/S has to appear as the expo- 
nent of a dimensionless quantity (l/L) in order to keep 
the dimension of Sp{l) unchanged. In all experiments 
and simulations of the 3 — D NS equation L is finite. In- 
fact a grand challenge for analytic theories of turbulence 
is to show that finite limit for C.p exist for I ^ L ^ oo. 
Such a scheme has been successfully carried out ||^ for 
the passive scalar field advected by a random velocity 
field in three dimensions (the Kraichnan model [Q). It 
has been shown that the anomalous scaling exponents 
of the passive scalar structure functions are independent 
of L (as L — > cxd), but the amplitudes do depend on L. 
Given that so much effort have been made to understand 
Sp{l,0), not much is known about Sp{l,t) even for small 
integer values of p (of course Sp{l, t) is more complicated 



than Sp{l,0)). Literature on Sp{l,t) in NS turbulence or 
related Burgers turbulence is rather sparse 0. Dynamic 
renormalisation group (DRG) calculations pi, one loop 
self-consistent calculations with the randomly forced 
3 — D NS equation suggest a long time dynamic expo- 
nent z = 2/3. In Ref [Q| large wavenumber limit of the 
velocity correlation function was also explored assuming 
dynamic scaling hypothesis to be valid. In experiments 
with a mean flow, because of large scale background ve- 
locity, one expects to measure z — 1. But even with zero 
mean velocity it has been shown , in the context of an 
1 — D Burgers equation, how z — I could arise. 

In this work we focus on the simplest two point non- 
equal-time velocity correlation function in the wave- 
number (fc) space. We show that in 3 — D fluid tur- 
bulence, in the large L limit, the real part of the 
non-equal time velocity correlation function Cij (k, t) = 
R{vi{'k, t)vj{~'k, 0)), for the inertial scales (/ = k~^), has 
oscillatory behaviour within a decaying envelope. From 
the incompressibility (V.v = 0) and isotropy assump- 



tions it follows Cij (k, t) 



c(fc,i)Fy(k). Here P^j{k) 
kikj /k"^ is the tranverse projector. 



An attempt to calculate the c{k, t) has been carried out 
by L'vov et.al. in the context of a turbulent flow with 
a mean velocity field Vq. But they treated Navier-Stokes 
equation at a linear level. They had compensated for the 
non-linear term to some extent by using a wave number 
dependent renormalised viscosity instead of the bare vis- 
cosity. They predict an oscillatory behaviour for c(fc,i), 
but in the absence of the mean velocity (i.e., Vo — 0) the 
oscillation vanish (i.e., it is a purely kinematic effect). 
But our data from a numerical simulation of the 3 — D 
NS equation (with large but finite L and zero mean veloc- 
ity) clearly reveals presence of oscillations. The data (see 
Fig. 4, 5) looks like the displacement of an under-damped 
harmonic oscillator. Also simulation of the REWA (re- 
duced wave vector set approximation) model by Eggers 
shows a non exponential decay at short times and a 
negative minima. 

Incompressible NS equation, forced randomly with a 
scale dependent variance, has been shown [ p^jjl^ to be 
a good model for fluid turbulence as far as multiscaling 
properties are concerned. But there exist many unre- 
solved theoretical problems with the analytic calculations 
with this model. Our calculation is based on a variant of 
this model, where instead of a singular forcing spectrum 
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(which goes as k~^) we use a spectrum which peaks at a 
smah but finite L~^, goes to zero a.t k — Q and behaves 
as for k ^ L~^. The equation of motion for the 
velocity field fourier component Vi(k) is 

w,(k) + i^ok^v^{k) = - iXM^ji{k) ^ Vj{q)vi{k - q) 

q 

+ MKt). (1) 

The random force fi(k,t) is a gaussian, white noise with 
the variancne 

(/,(k,i)/,(k',t')) = ^^^^^n,imk+k')6it~t') 

(2) 

Here Miji = [fcjP,;;(k) + fciPij(k)]/2 and A is an artificial 
coupling constant which will be set to 1 later. Henceforth 
we will denote the variance of the force 2Dok/{k'^ + L^^)'^ 
by D{k). 

In the theories with singular forcing spectrum (or 
equivalently infinite integral scale), infrared divergences 
appear if one tries to calculate effective viscosity per- 
turbatively. Also the scheme cannot handle the so called 
sweeping effect i.e., the interaction of the bigger eddies (of 
size q~^) with the eddy of size k~^ (when {q < k)). Phys- 
ically the bigger eddies just advect the smaller eddies 
without distorting them much, so such divergences are 
basically defect of such a pcrturbative scheme. But One 
should remember that this eddy picture is quite heuristic 
in nature because velocity fourier modes v{k) are global 
features of the velocity field where as the eddies are spa- 
tially correlated patches in the velocity field and hence 
local in nature. One systematic way to get rid of sweeping 
divergences in the equal time velocity structure functions, 
is to go to the lagrangian frame. Another way is to do an 
RG calculation which excludes the effect of the {q < k) 
modes on the k mode, so both the infrared divergence 
and sweeping effect are eliminated. In these calculations 
because of universality reasons one is mainly interested 
in the zero frequency limit of the effective viscosity i.e., 
5i'{k,uj — > 0) and assumes that for all frequencies to, the 
effective propagator G{k,uj) ^ {~iuj + k'^5v{k,uj 0))^^ 
i.e., remains a Lorentzian in oj. This approximation 
works well for long time properties. But here since we 
are interested in c(fc, t) for all t (including the short time 
behaviour) we need the correct behaviour of G{k,uj) for 
all Lus'. 

Our procedure to calculate c{k,Lu) is a mixture of a 
self-consistent and a perturbative scheme. We show that 
if we assume a large but finite integral scale L, due to 
nonlinear interaction among modes, the effective viscos- 
ity is complex. It has the regular renormalized real part 
and a k,uj dependent imaginary part as well. The os- 
cillation in c{k,t) arises because of this imaginary part. 



We use an one loop perturbation theory to determine the 
complex viscosity. The calculation is similar to the stan- 
dard RG procedure for evaluating zero frequency viscos- 
ity. But unlike in the RG procedure, where only modes 
greater than the external wave-vector k are integrated 
out, we integrate over all q modes, including the range 
[0, k]. In the zero frequency limit {co = 0) our calculation 
is self-consistent (at one loop level), but for finite ui it is 
a perturbative calculation. Since our forcing spectrum is 
not singular, there is no infrared divergence in our inte- 
grals. 

Treating the nonlinear term perturbatively |§] the ef- 
fective response function G{k,uj) can be calculated as 

= —ioj + v^k^ + 5v{k^ oj). In the small k and vq ^ Q 
limit, v^k^ is negligible compared to 5v{k,uj). Hence 

= -i[uj - k'^I{S)] + k'^R{S), where we have denoted 
5v{k,uj) by and the real, imaginary parts by R,I. 

Using c{k,Lo) = D{k)\G{k,uj)\'^ we get 

c(fc,w) = 2 (3) 

[Lo~k-^I{S)Y + k^R{Sf 

c{k,t) is the inverse fourier transform of c(fc,a;). So 
our task is to calculate S. Following Ref. |^ 

-PSPl,{k) = [iXfMirnnO^) J ^^^M„,,(k-q) X 

P,„M)D{q)\G{q,cj')\'Gi\k~q\,u;-u;') 

(4) 

Multiplying both sides by Pji (k) and contracting over l,j 
we get 

" 2fc2 J "(2^ " q)P^mic^)Diq) x 

|G(g,a;')pG(|k-q|,c.-c.') 
Integrating the r.h.s. over u' gives 

/ ,„ ,„ o(k, k - q, q)- -— (5) 

J (27r)3 ^ ^'^^2zyig^ -zcj + i/i(g^ + |k-q|^) ^ ' 

k^b{k, k — q, q) is obtained by contracting Mjmn, Mnij 
and Pirn ■ The expression for b is b{k, k ~ q, q) = 
^^-^ (xy + z^) llj] . The trio (k, k — q, q) form a triangle 
and x,y,z are the direction cosines of the angles oppo- 
site to k, k — q, and q respectively. We have calculated 
this integral numerically. We have used G{k,u!)~^ = 
—iLu + vik^ anticipating a renormalisation of the viscos- 
ity in the w limit. Here vi could be a function of 
L, Da which we determine later. But as far as the fc, lo 
dependence of the integral is concerned, for fixed L and 
Dq, we can treat vi as a constant. From Eq.5 note that 
I{S) — > 0, as w ^ 0. We check numerically that at w = 0, 
the integral in Eq.5 scales as fc~^^. Also by expanding 
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FIG. 1. viiL) versus L. It shows convergence of vi{L) 
as L — > oo. In our numerical evaluation of the integral in 
Eq.5 we choose Dq — 0.01 and first evaluate the integral for 

= in order to get Fig.l. From this and Eq.6 we determine 
vi ~ 0.06. Then we choose a large enough value for L{= 100) 
and evaluate the integral for nonzero values of u to get S. 

the integral in small k and retaining only leading powers 
in k (as is done in a DRG procedure ^) one can infer 
this. If self-consistency has to be achieved at w = 0, then 
on the l.h.s. of Eq.5, i?[S'(fc,0)] = j/ifc^-^. This fixes 
z = 2/3. Numerical evaluation of the integral with a fi- 
nite L corroborates the approximate analytic prediction 
for small k because L~'^ « k. Now we determine vi in 
a self consistent way Q . 

At w = 0, self consistency of Eq.5 requires 

(6) 

Where g{L) is the L dependent integral on the r.h.s. 
of Eq.5. In Fig.l we show that g{L) converges to a finite 
value as L grows large. Using g{L) ~ .95 from Fig.l, 

1 /3 

we get viiL) ~ 0.28Do • While evaluating the integral 
for inertial range k modes we used G^^ = —iuj + vik'^^^ 
(neglecting the v^k^ term) for all the q and k — q modes 
in the integrand, though the integral runs over both in- 
ertial and dissipation modes. But this is approximately 
correct because modes very far from the k mode do not 
contribute (locality in k space). 

Returning to Eq.5, the dominant pole of c{k,uj) will 
decide the oscillation and the decay of c{k,t). If the 



FIG. 2. S versus lj for two different k values. Real and 
imaginary parts of S are indicated by R[S] and I[S]. 

pole lies at w = cji + iuj2 then lui ~ k'^I{S[uJi]) and 
= -k^R{S[Lj2]). From the shapes of R, I{S) in Fig.2 
and the fact that R, I{S) are even, odd functions of ui, we 
can infer that there will be two solutions to these tran- 
scendental equations. They are of the form ±wi + 1^2, 
when UJ2 is negative. This ensures the causality of the 
effective response function G. The approximate data col- 
lapse for the inertial range k modes (?7^^ > fc > L^^) in 
our Fig. 3 implies that the oscillation period scales with 

fc-2/3. 

We performed a DNS of the 3 — D NS equation and 
calculated c{k,t) versus t for various v(k) modes. The 
data (see Fig. 4, 5) clearly shows oscillation in time. In 
our simulation forcing was present only at large length 
scales (i.e., in the fc-space the v(k) modes in the smallest 
two shells were forced). Our pseudo-spectral scheme for 
the DNS is same as in . We used a 32"^ grid with pe- 
riodic boundary condition. We obtained a short inertial 
range (shown in Fig. 6) with Re\ ~ 22. We have obtained 
very long time series (> 68T{L)) for averaging c{k,t). A 
preliminary test run |T7| on a 64'^ grid also confirms ex- 
istence of such oscillations. 

Our simulation data in Fig. 4 indicates that the oscil- 
lation time period approximately scales as k~^^^ in the 
inertial range. But Fig. 5 shows that modes closed to 
the forced modes and the dissipative modes differ widely 
from the inertial ones in periodicity and decay. We have 
averaged the data for > 68T(L) (when T{L) is the large 
eddy turnover time). This is sufhcienty long averaging 
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FIG. 3. c{k,t) versus the scaled time t/T{k) (with 
T{k) = k~'^^^), for two different k values. 



FIG. 5. Cij{k,t)/Cij{k,0) versus t/T{k) plotted for three 
different k vectors which lie close to the forced range (cir- 
cle), in the inertial range (square) and close to the dissipation 
range (triangle). They show wide separation in periodicity 
and decay. 




15 20 25 30 35 40 
t/T(k) 

FIG. 4. Cij(k,t)/C'ij(k,0) plotted against scaled time 
t/T{k) (with eddy turnover time T{k) = Ae-'^^^k'^^^). Plots 
for different, independent k vectors (in the inertial range) and 
different i,j are indicated by k = {kx,ky,k^) : Approx- 
imate data collapse can be seen. Data point shown at the 
longest time separation and zero time separation are aver- 
aged over ~ 68 and ~ 84 Large-eddy turnover times (T(L)) 
respectively. 
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Energy spectrum E(k) 
■ slope = -5/3 
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FIG. 6. The energy spectrum InE(k) versus Ink. The po- 
sition of the different chosen k values, for which c(k, t) has 
been plotted in Fig. 4 and 5, are indicated at the top of the 
figure. 
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time for an equal time, shell averaged, correlation func- 
tions to converge; but here for c(k,t) since (a) single k 
mode is involved, (b) long time history (few T(L)'s) of the 
mode is important, it requires a longer averaging time. 
In a simulation with small but finite the vik^^^ term 
will loose its dominance over fofc^ as k increases towards 
the dissipation range. That explains the strong damping 
seen in the k = (0, 5, —6) mode in Fig. 5. 

In the simulation we calculate T(L) using the for- 
mula T{L) = Lbox/vrms, whcre Lbox is the simula- 
tion box size. From dimensional analysis arguments 
1^ the scale dependent eddy-turnover time T{k) = 
Ae~^^^k~^^^ , when A is a constant of 0{1). Equating 
T{L) = yle-i/3(27r/L)-2/3 ^e get the prefactor Ae^i/^ 
and hence can calculate T(fc). Fig. 4 shows that the time 
period of oscillation Afc ~ 25T{k). This ratio may be non- 
universal. To get a clue let us look at our randomly forced 
model, where I{S) is an explicit function of Dq, i^i and 

1 /3 

vi oc Dq (see Eq.5 and 6). Hence Afc depends on Dq in a 
complicated way. We explore the dependence of T{k) on 
Dq below. In this model T(k) cannot be determined in a 
simple way because here all the k shells are being forced 
and hence the energy flux is not a constant but increses 
logarithmically with k The scale dependent energy 
fluxn(fc) = Diq)d^q/{2Trf ^ £^[\n{a^ + e)-k''{a^ + 
k'^)~^ — 21n(a)], when a = L^^. Again using dimen- 
sional analysis we get T{k) = A.Il{k)^^/^k^^/^ (hence 
T{k) cx D~^^^). Evaluating U{k)-^/^ for Dq = 0.01 
and L — 100 (which we used for our theoretical graphs 
in Fig.2,3) gives U{0.1)-^^^ = 6.17, n{1.0)-^/^ = 8.1. 
AUthough this k dependence is weak, but the perfect 
data collapse with respect to tk^^^ in Fig. 3 implies that 
the above dimensional analysis estimate is not accurate. 
In Fig. 3 the oscillation period ^ G/c^^/'^ and hence 
Afe ~ T{k) (neglecting the constant A of 0(1)). 

The imaginary part of the viscosity, generated by an 
one loop perturbation theory, cannot be interpreted as 
a background velocity which is slowly varying in time. 
This mis-interpretation may be provoked by the fact 
that if we had a mean background flow Vq in the prob- 
lem, then the nonlinear term would generate an extra 
term ?(Vo.k)v(k, in the NS equation. With this ex- 
tra linear term the bare propagator will be G{k,uj) ~ 



(- 



iVo.k + i^ok 



Il0f|. But the k,uj dependent 



complex viscosity, which we find in our theory, when 
transformed to (x, t) space, gives an additional memory 
term J dx'dt'K(|x — x'|,< — <')v(x',<') in the equation of 
motion (e.o.m.). But it is true that the K{k,t) field oscil- 
lates at a slower time scale than that of the v(k, t) mode 
itself (as I5v{k,uj) has a peak at a lower lo than that of 
c(k, Lo)). Also a significant contribution to the integral for 
5v{k, oj) comes from the q < k modes. So it does resemble 
sweeping by larger (and hence slower) eddies to some ex- 



tent. But rigorously convection by a background velocity 
Vo(x, t) should look like Vo(x, i).Vv(x, t) (which is local 
in (x, i)). We note that in the turbulence context, com- 
plex effective viscosity has been proposed before by J.K. 
Bhattacharjee in Q. The author had assumed dynamic 
scaling hypothesis (as in dynamic critical phenomena) to 
be valid for fluid turbulence and had predicted a form for 
the effective viscosity in the high frequency limit. Then 
an interpolation scheme had been used to connect the 
two limits (small uj and large oj) of the effective viscosity. 

In conclusion, we have given numerical evidence and 
an approximate theory for the novel oscillatory decay of 
the two-point, temporal correlation function in 3 — D 
fluid turbulence. This behaviour is similar to viscoelastic 
effect seen in complex fluids. 
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